\documentclass[a4paper, 14pt]{article}
\usepackage{CJKutf8}
\usepackage{graphicx}
\usepackage{listings}
\lstset{breaklines}
\lstset{extendedchars=false}
\lstset{language=C}
\usepackage{graphicx}
\usepackage{pythonhighlight}
\usepackage{amsmath,amsfonts,amsthm,amssymb}
\theoremstyle{definition}
\newtheorem{thm}{Theorem}
\newtheorem{exmp}{Example}
\newtheorem{defn}{Definition}
\newtheorem{lema}{Lemma}
\newtheorem{prop}{Proposition}
\newtheorem{coro}{Corollary}


\renewcommand{\baselinestretch}{1.25}


%------------setlength----------------%
\setlength{\textwidth}{162mm}
%\setlength{\textheight}{190mm}
\setlength{\textheight}{231mm}
\setlength{\headheight}{-0.1cm}
\setlength{\topmargin}{-0.1cm}
\setlength{\oddsidemargin}{-0cm}
\setlength{\evensidemargin}{-0cm}

\setlength{\parskip}{1mm}
\setlength{\unitlength}{1mm}
\setlength{\parindent}{2em}

\title{用最小二乘法求解蒸汽和液体的平衡问题}

\author{李之琪\\数学科学学院，数学与应用数学\\学号：3180103041}


\begin{document}
\begin{CJK*}{UTF8}{gbsn}
\maketitle
\section{问题}

蒸馏工艺中，我们通过混合物各组分的沸点不同，用不同温度的蒸发和冷凝分离组分。其中，描述二元气-液平衡状态的数学模型是：
$$
x_{i} \nu_{i} p_{i}^{\mathrm{sat}}=y_{i} p, i=1,2
$$
这里，$p$是总压，$x_{i}$和$y_{i}$分别是组分$i$的液相和气相摩尔分数，$p_{i}^{\mathrm{sat}}$和$\nu_{i}$分别是组分$i$在平衡温度$T$时的饱和蒸汽压和活度系数，$i=1,2$。整理后，得到：
$$
p=x_{1} \nu_{1} p_{1}^{\mathrm{sat}}+\left(1-x_{1}\right) \nu_{2} p_{2}^{\mathrm{sat}}
$$
其中，$p_{i}^{\mathrm{sat}}$满足Antoine方程：
$$
\log_{10} p_{i}^{\mathrm{sat}}=c_{1 i}-\frac{c_{2 i}}{T+c_{3 i}}
$$
$c_{1 i}$，$c_{2 i}$和$c_{3 i}$是已知系数，$i=1,2$。二元的$\nu_{i}$满足Van Laar模型：
$$
\begin{array}{l}
\ln \nu_{1}=a_{1}\left(\frac{a_{2} x_{2}}{a_{1} x_{1}+a_{2} x_{2}}\right)^{2} \\
\ln \nu_{2}=a_{2}\left(\frac{a_{1} x_{1}}{a_{a} x_{1}+a_{2} x_{2}}\right)^{2}
\end{array}
$$
其中$a_1$，$a_2$是待定系数。
给定一组水和(1,4)二氧杂环乙烷在$T=20$时的气液平衡实验数据$(x_{1 i},p_i)$，$i=1,2,\cdots,11$。根据上述数据和模型，确定公式中的$a_1$，$a_2$，从而给出确定的水和(1,4)二氧杂环乙烷在$T=20$时的气液平衡公式。

\section{模型建立}
我们采用最小二乘模型确定公式中的待定系数$a_1$和$a_2$。记
$$
p\left(\textbf{x}, a_{1}, a_{2}\right)=x_{1} \exp \left(a_{1}\left(\frac{a_{2} x_{2}}{a_{1} x_{1}+a_{2} x_{2}}\right)^{2}\right) p_{1}^{\mathrm{sat}}+x_2 \exp \left(a_{2}\left(\frac{a_{1} x_{1}}{a_{1} x_{1}+a_{2} x_{2}}\right)^{2}\right) p_{2}^{\mathrm{sat}}
$$
这里$\textbf{x}=(x_1,x_2)$。由于$x_1+x_2=1$，上式本质上只有一个自由变量，记为$p\left(x_1,a_1,a_2\right)$。则残量:
$$
r_j(a_1,a_2) = p\left(x_{1 j},a_1,a_2\right) - p_j
$$
最小二乘模型的目标函数为：
$$
f(a_1,a_2)=\frac{1}{2}\sum_{j=1}^{m} r^2_j(a_1,a_2) =
\frac{1}{2}\sum_{j=1}^{m}[p\left(x_{1 j},a_1,a_2\right) - p_j]^2
$$
这里$m=11$。我们的目标是寻找$f(a_1,a_2)$的最小值点。这是一个非线性最小二乘问题，记
$$
\begin{array}{l}
r(\textbf{a})=(r_1(\textbf{a}),r_2(\textbf{a}),\cdots,r_m(\textbf{a}))^T \\
J(\textbf{a})=(\frac{\partial r_j}{\partial a_i} )_{j=1,2,\cdots,m,i=1,2}
\end{array}
$$
其中$\textbf{a}=(a_1,a_2)$。舍去高阶项，可以将上述非线性最小二乘问题转化为以下一系列信任域问题的求解：
$$
\min_p m_k(p)=\min_p \frac{1}{2}\left \|J_kp+r_k  \right \|^2_2, s.t.
\left \| p \right \| \le \triangle k 
$$
这里$\triangle k$是信任域半径。我们将采用LM算法求解该信任域问题，从而得到最小二乘模型的解。

\section{算法构建}
将问题化为信任域问题的求解后，我们采用LM算法进行求解。LM算法的思路较简单：在第$k$步，对某个$\lambda_k$，求解处理后的方程
$$
(J^T_kJ_k+\lambda_k I)p_k = -J^T_kr_k
$$
根据$p_k$，计算$\rho_k=\frac{\triangle f_k}{\triangle m_k}$，用来衡量$p_k$的效果。若$\rho_k$太小，说明此步所求的$p_k$对模型下降效果较差，此时我们放大$\lambda_k$，这实质上是缩小了信任域半径$\triangle k$，使得下一步在更小更精细的范围内求解；若$\rho_k$接近于1，说明此步所求的$p_k$对模型下降效果很好，下降速度被信任域半径限制住了，此时我们缩小$\lambda_k$从而放大信任域半径$\triangle k$，使得下一步在更大范围内求解。在$\lambda_k$的调整结束后，我们根据下降效果的好坏决定是否前进到$\textbf{a}_{k+1}=\textbf{a}_k + p_k$。\\
如此重复以上步骤，直至目标函数$f(\textbf{a})$充分小或$k$达到指定的循环次数。最终得到的$\textbf{a}=(a_1,a_2)$即为最小二乘模型的求解结果。\\
我们使用2011标准的C++编写程序，下面将结合部分代码对LM算法进行进一步的说明。首先对项目文件的结构和功能做简单的说明：\\\\
\textbf{Project}\\
\textbf{——LSM.h}:包含所有类和函数定义，是程序的主体部分。具体又包括以下几个类：\\
\textbf{————class function}：抽象的模板类，用于表达一般的从$\mathbb{R}^n$到$\mathbb{R}$的函数，其主要成员函数如求值、求导数用纯虚函数提供接口，供子类实现。\\
\textbf{————class GLB\_function}：由function$<$2$>$派生得到的子类，提供函数$p(x,a_1,a_2)$及其导数的实现。\\
\textbf{————class LSM}：生成最小二乘方法的模板类，将类GLB\_function和提供的多组数据作为其成员变量。提供LM算法用于求解最小二乘问题。\\
\textbf{——main.cpp}：从文件Inputfile中读取初始参数、期望精度等变量，把这些变量作为函数参数调用LM算法求解问题，输出结果并得到文件result.m，用matlab运行该文件得到可视化结果。\\
\textbf{——Antoine-coefficients}：存储Antoine方程中需要用到的系数，在class GLB\_function中被读取。\\
\textbf{——data}：存储用于构建最小二乘模型的数据，在class LSM中被读取。\\
\textbf{——Inputfile}：存储初始参数、期望精度、初始$\lambda$，在main.cpp中被读取。\\
\textbf{——Makefile}：通过命令行指令"make"得到可执行文件main；"make test"得到可执行文件test。\\
\newpage
\noindent 下面展示用C++实现的LM算法：
\begin{lstlisting}
template<int n_coe, int n_data>
double* LSM<n_coe,n_data>::LM(double* _coes, double tol , double lamda){
  double f=this->f(this->r(_coes));
  int k = 0;
  while (f > tol && k < 1000){
	Matrix<double,n_data,n_coe> J = this->J(_coes);
	Matrix<double,n_coe,n_coe> A = J.transpose()*J+lamda*MatrixXd::Identity(n_coe, n_coe);
	Vector<double,n_coe> df = this->df(_coes);
	Vector<double,n_coe> p = A.colPivHouseholderQr().solve(-df);
	double rho = this->rho(_coes,p,lamda);
	if (rho < 0.25) lamda = 4 * lamda;
	if (rho > 0.75) lamda = 0.5 * lamda;
	if (rho > 0.1){
	  for (int i = 0 ; i < n_coe ; i++)
		_coes[i]+=p(i);
	  f = this->f(this->r(_coes));
	}
	k++;
  }
  return _coes;
}
\end{lstlisting}
值得一提的是，由于本问题需要求解的线性方程组的矩阵规模仅为$2\times2$,求解采用了简单的基于Householder变换的QR分解。若矩阵规模较大，这里应该使用Givens变换以降低运算量。

\section{结果展示}
通过命令行指令"make"得到可执行文件main，运行main即可得到计算结果。下面展示计算结果。输入文件Inputfile如下：
\begin{table}[!htp]
	\centering
	\begin{tabular}{|c|c|c|c|}
		\hline
		initial of $a_1$ & initial of $a_2$  & tol & lambda  \\
		\hline	
		1     & 1            &0.1 & 5000           \\
		\hline
		
	\end{tabular}
\end{table}\\
结果为:\\
a1 = 1.95842 , a2 = 1.68919\\
即水和(1,4)二氧杂环乙烷在$T=20$时的气液平衡公式为：
$$
p\left(\textbf{x}\right)=x_{1} \exp \left(1.95842\left(\frac{1.68919 x_{2}}{1.95842 x_{1}+1.68919 x_{2}}\right)^{2}\right) p_{1}^{\mathrm{sat}}+x_2\exp \left(1.68919\left(\frac{1.95842 x_{1}}{1.95842 x_{1}+1.68919 x_{2}}\right)^{2}\right) p_{2}^{\mathrm{sat}}
$$
这里$\textbf{x}=(x_1,x_2)$，$p_{1}^{\mathrm{sat}}$，$p_{2}^{\mathrm{sat}}$可由Antoine方程得到。
\newpage
\noindent 我们来观察公式与实际数据的偏差。下面列出不同$x_1$时用以上公式算出的结果和实际的结果的比较。以下结果可以通过命令行"make test"产生可执行文件test并运行之得到。
\begin{table}[!htp]
	\centering
	\begin{tabular}{|c|c|}
\hline
result by formula & result from datas\\
\hline
28.8241        &      28.1\\
\hline
34.6443        &      34.4\\
\hline
36.453        &      36.7\\
\hline
36.8673        &      36.9\\
\hline
36.874        &      36.8\\
\hline
36.7498        &      36.7\\
\hline
36.3904        &      36.5\\
\hline
35.3848        &      35.4\\
\hline
32.9478        &      32.9\\
\hline
27.73        &      27.7\\
\hline
17.4733        &      17.5\\
\hline
	\end{tabular}	
\end{table}\\
可视化结果如下：(用matlab运行result.m)
\begin{figure}[!htp]   
	\centering
	\includegraphics[width=10cm]{1.eps}
\end{figure}\\
其中，蓝线表示公式的结果，红色星号表示用于计算的数据。可以看出，该公式效果整体效果较好；更具体而言，在$x$较小时($x \le 0.3$)，误差要相对大一些；在$x$较大时，误差较小。
\section{总结}
在本次项目作业中，我们针对蒸汽和液体的平衡问题构建了最小二乘模型，并使用C++利用LM算法求解该模型，得到了指定温度下的气液平衡公式。最后，我们将该公式的计算结果与实际数据相比较，发现公式的拟合效果整体上较好，误差不大。
\end{CJK*}
\end{document}